/********************************************************************
*  sample.cu
*  This is a example of the CUDA program.
*********************************************************************/
#include <terrainCuda.h>



/************************************************************************/
/* Init CUDA                                                            */
/************************************************************************/
#if __DEVICE_EMULATION__

bool InitCUDATerrain(void){return true;}

#else
bool InitCUDATerrain()
{
	int count = 0;
	int i = 0;

	cudaGetDeviceCount(&count);


	for(i = 0; i < count; i++) {
		cudaDeviceProp prop;
		if(cudaGetDeviceProperties(&prop, i) == cudaSuccess) {
			if(prop.major >= 1) {
				break;
			}
		}
	}

	cudaSetDevice(i);


	return true;
}

#endif
/************************************************************************/
/* Functions  - Perlin noise                                                           */
/************************************************************************/



__device__ void normalize (float3 * a)
{
	float l = sqrt(a->x*a->x + a->y*a->y + a->z*a->z);
	a->x/=l;
	a->y/=l;
	a->z/=l;
}

__device__ float Noise(int x, int y)
{
	int n = x + y * 57;
	n = (n<<13) ^ n;
	return  (1.0 - ( (n * (n * n * 15731 + 789221) + 1376312589) & 0x7fffffff) / 1073741824.0);    
}

__device__ float SmoothedNoise1(float x, float y)
{
	float corners = ( Noise(x-1, y-1)+Noise(x+1, y-1)+Noise(x-1, y+1)+Noise(x+1, y+1) ) / 16;
	float sides   = ( Noise(x-1, y)  +Noise(x+1, y)  +Noise(x, y-1)  +Noise(x, y+1) ) /  8;
	float center  =  Noise(x, y) / 4;
	return corners + sides + center;
}

__device__ float Interpolate(float a, float b, float x)
{
	float ft = x * 3.1415927;
	float f = (1 - cos(ft)) * 0.5;

	return  a*(1-f) + b*f;

}




__device__  float InterpolatedNoise1(float x, float y)
{

	int integer_X    = int(x);
	float fractional_X = x - integer_X;

	int integer_Y    = int(y);
	float fractional_Y = y - integer_Y;

	float v1 = SmoothedNoise1(integer_X,     integer_Y);
	float v2 = SmoothedNoise1(integer_X + 1, integer_Y);
	float v3 = SmoothedNoise1(integer_X,     integer_Y + 1);
	float v4 = SmoothedNoise1(integer_X + 1, integer_Y + 1);

	float i1 = Interpolate(v1 , v2 , fractional_X);
	float i2 = Interpolate(v3 , v4 , fractional_X);

	return Interpolate(i1 , i2 , fractional_Y);
}



__device__  float PerlinNoise_2D(float x, float y, float persistence, int octaves)
{

	float total = 0;
	float p = persistence;
	int n = octaves;

	for(int i=0; i<n; i++)
	{

		float frequency = pow(2.0f,(float)i);
		float amplitude = p*i;

		total = total + InterpolatedNoise1(x * frequency, y * frequency) * amplitude;
	}



	return total;
}



__global__ void generate(SimpleVertex *a, int N, float scale)
{
	int idx = blockIdx.x * blockDim.x + threadIdx.x;

	int X= idx%N;
	int Y= (int)((idx-X)/N);



	//Calculate position
	a[idx].Pos.x=X*scale;
	a[idx].Pos.y=PerlinNoise_2D(((float)X)/20,((float)Y)/20, 1, 4)*scale;
	a[idx].Pos.z=Y*scale;

	//Calculate texture coordinates
	a[idx].Tex.x=X* 0.1; 
	a[idx].Tex.y=Y* 0.1; 

}


///Normal generation

__device__ float3 cross(float3 a, float3 b)
{
	float3 c=make_float3(a.y*b.z-a.z*b.y,
		a.z*b.x-a.x*b.z,
		a.x*b.y-a.y*b.x);

	return c;
}


__global__ void generateNormals(SimpleVertex *t, int N)
{
	int idx = blockIdx.x * blockDim.x + threadIdx.x;

	int X= idx%N;
	int Y= (int)((idx-X)/N);

	float3 normal;

	if(X>0&&Y>0&&X<N-1&&Y<N-1)
	{
		float3 toA=make_float3(t[idx-N].Pos.x-t[idx].Pos.x, t[idx-N].Pos.y-t[idx].Pos.y, t[idx-N].Pos.z-t[idx].Pos.z);
		float3 toB=make_float3(t[idx+1].Pos.x-t[idx].Pos.x, t[idx+1].Pos.y-t[idx].Pos.y, t[idx+1].Pos.z-t[idx].Pos.z);
		float3 toC=make_float3(t[idx+N].Pos.x-t[idx].Pos.x, t[idx+N].Pos.y-t[idx].Pos.y, t[idx+N].Pos.z-t[idx].Pos.z);
		float3 toD=make_float3(t[idx-1].Pos.x-t[idx].Pos.x, t[idx-1].Pos.y-t[idx].Pos.y, t[idx-1].Pos.z-t[idx].Pos.z);

		float3 AB=cross(toB,toA);
		float3 BC=cross(toC,toB);
		float3 CD=cross(toD,toC);
		float3 DA=cross(toA,toD);

		normal=make_float3((AB.x+BC.x+CD.x+DA.x)/4,
			(AB.y+BC.y+CD.y+DA.y)/4,
			(AB.z+BC.z+CD.z+DA.z)/4);


		float l=sqrt(normal.x*normal.x+normal.y*normal.y+normal.z*normal.z);

		t[idx].Nor.x=normal.x/l;
		t[idx].Nor.y=normal.y/l;
		t[idx].Nor.z=normal.z/l;

		normalize(&normal);
		float3 tangent= cross(make_float3(0,normal.y,normal.z), make_float3(1,0,0));
		normalize(&tangent);
		float3 binormal= cross(tangent,make_float3(normal.x,normal.y,0));

		t[idx].Tan.x=tangent.x;
		t[idx].Tan.y=tangent.y;
		t[idx].Tan.z=tangent.z;

		t[idx].Bin.x=binormal.x;
		t[idx].Bin.y=binormal.y;
		t[idx].Bin.z=binormal.z;

	}

	else
	{
		t[idx].Nor.x=0.0f;
		t[idx].Nor.y=1.0f;
		t[idx].Nor.z=0.0f;

		normal=make_float3(0,1,0);



		normalize(&normal);
		float3 tangent= cross(make_float3(0,normal.y,normal.z), make_float3(1,0,0));
		normalize(&tangent);
		float3 binormal= cross(tangent,make_float3(normal.x,normal.y,0));
		normalize(&binormal);

		t[idx].Tan.x=tangent.x;
		t[idx].Tan.y=tangent.y;
		t[idx].Tan.z=tangent.z;

		t[idx].Bin.x=binormal.x;
		t[idx].Bin.y=binormal.y;
		t[idx].Bin.z=binormal.z;
	}



}



/************************************************************************/
/* Main                                                           */
/************************************************************************/
SimpleVertex* generateTerrain(int s, float scale,  ofstream &out)
{

	if(!InitCUDATerrain()) {
		return 0;
	}

	SimpleVertex *a_h, *a_d;  // Pointer to host & device arrays
	const int N = s*s;  // Number of elements in arrays
	size_t size = N * sizeof(SimpleVertex);

	a_h = (SimpleVertex *)malloc(size);        // Allocate array on host
	cudaMalloc((void **) &a_d, size);   // Allocate array on device

	// Initialize host array and copy it to CUDA device

	cudaMemcpy(a_d, a_h, size, cudaMemcpyHostToDevice);

	// Do calculation on device:
	int block_size = 16;
	int n_blocks = N/block_size + (N%block_size == 0 ? 0:1);
	generate <<< n_blocks, block_size >>> (a_d, s, scale );
	cudaThreadSynchronize();
	generateNormals <<< n_blocks, block_size >>> (a_d, s);
	cudaThreadSynchronize();

	// Retrieve result from device and store it in host array
	cudaMemcpy(a_h, a_d, sizeof(SimpleVertex)*N, cudaMemcpyDeviceToHost);

	// Print results
	for (int i=0; i<N; i++) out<< i<<" "<< a_h[i].Tex.x<<" "<<a_h[i].Tex.y<<" "<<n_blocks<<endl;;

	// Cleanup
	cudaFree(a_d);

	return a_h;
}


SimpleVertex* generateTerrain(int s, float scale )
{

	if(!InitCUDATerrain()) {
		return 0;
	}

	SimpleVertex *a_h, *a_d;  // Pointer to host & device arrays
	const int N = s*s;  // Number of elements in arrays
	size_t size = N * sizeof(SimpleVertex);

	a_h = (SimpleVertex *)malloc(size);        // Allocate array on host
	cudaMalloc((void **) &a_d, size);   // Allocate array on device

	// Initialize host array and copy it to CUDA device

	cudaMemcpy(a_d, a_h, size, cudaMemcpyHostToDevice);

	// Do calculation on device:
	int block_size = 16;
	int n_blocks = N/block_size + (N%block_size == 0 ? 0:1);
	generate <<< n_blocks, block_size >>> (a_d, s, scale );
	cudaThreadSynchronize();
	generateNormals <<< n_blocks, block_size >>> (a_d, s);
	cudaThreadSynchronize();

	// Retrieve result from device and store it in host array
	cudaMemcpy(a_h, a_d, sizeof(SimpleVertex)*N, cudaMemcpyDeviceToHost);



	// Cleanup
	cudaFree(a_d);

	return a_h;
}

